use 5.10.0;
use strict;

package Graphics::PlanarVector::PathSegment;
use List::Util qw/min max/;
use Graphics::PlanarVector::Vector qw/vector veq vdiv/;
Math::Util::Range->import('range');
for (qw/poly roots/) {
	Math::Util::Polynomial->import($_)
}

# pt: points as in bezier curve
# poly: polynomial form in terms of t
use fields qw/pt poly uid bbox rev sect rem depth parent/;
use Carp;
use Exporter 'import';
our @EXPORT_OK = qw/path_segment/;

use overload (
	'""'	=> sub {
		my $s = shift;
		local $" = ', ';
		"seg($s->{uid})[@{$s->{pt}}] ".$s->{parent}[0]
	},
);

my $uid = 0;

sub path_segment {
	my $self = fields::new(__PACKAGE__);
	my @pts;
	while (@_) {
		my $arg = shift;
		push @pts, ref($arg) 	? $arg
					: vector($arg, shift)
	}
	$self->{pt} = \@pts;
	{
		lock $uid;
		$self->{uid} = ++$uid;
	}
	$self->{parent} = [$uid, 0, 1];
	return $self
}

sub dup {
	my $p = shift;
	my $r = path_segment(@{$p->{pt}});
	$r->{bbox} = $p->{bbox};
	$r->{parent} = [@{$p->{parent}}];
	$r
}

sub from_poly {
	my $p = shift;
	pop @$p while ($p->[-1]->len2 < 1e-10);
	return unless @$p >= 2;

	my @pts;
	if ($#$p == 3) {
		@pts =  ($p->[0], $p->[0] + $p->[1] / 3,
			$p->[0] + ($p->[1] * 2 + $p->[2]) / 3,
			$p->[0] + $p->[1] + $p->[2] + $p->[3]);
	}
	elsif ($#$p == 2) {
		@pts = ($p->[0],
			$p->[0] + $p->[1] / 2, 
			$p->[0] + $p->[1] + $p->[2]);
	}
	elsif ($#$p == 1) {
		@pts = ($p->[0], $p->[0] + $p->[1])
	}
	path_segment(@pts);
}

sub reverse {
	my $path = shift;
	$path->{rev} //= do {
		my $seg = path_segment(reverse @{$path->{pt}});
		@{$seg->{parent}} = @{$path->{parent}}[0,2,1];
		$seg;
	}
}

sub seg_equal {
	my ($this, $that) = @_;
	return unless order($this) == order($that);
	for (0 .. $#{$this->{pt}}) {
		return unless veq($this->{pt}[$_], $that->{pt}[$_])
	}
	return 1;
}

my %to_poly = (
	1 =>	sub {
		my ($p, $off) = (shift, shift // vector(0,0));
		my @s = map { $_ + $off } @{$p->{pt}};
		poly($s[0], $s[1] - $s[0])
	},
	2 =>	sub {
		my ($p, $off) = (shift, shift // vector(0,0));
		my @s = map { $_ + $off } @{$p->{pt}};
		poly($s[0], -2 * ($s[0] - $s[1]), $s[0] - 2 * $s[1] + $s[2])
	},
	3 =>	sub {
		my ($p, $off) = (shift, shift // vector(0,0));
		my @s = map { $_ + $off } @{$p->{pt}};
		poly(	$s[0],
			-3 * ($s[0] - $s[1]),
			3 * ($s[0] - 2 * $s[1] + $s[2]),
			$s[3] - $s[0] + 3 * ($s[1] - $s[2])
		)
	},
);

sub order($){ @{shift->{pt}} - 1 }

sub poly_form {
	my ($s, $off) = @_;
	if (!$off) {
		return $s->{poly} //= $to_poly{order($s)}->($s)
	}
	$to_poly{order($s)}->($s, $off)
}

sub param_to_point {
	my ($p, $t) = @_;
	my $pt = $p->{pt}[0];
	$t <= 0 ? $pt : $t >= 1 ? $p->{pt}[-1] : do {
		&poly_form($p, -$pt)->poly_eval($t) + $pt
	};
}

sub valid_param {
	my %seen;
	grep { !ref($_) &&  $_ >= 0 && $_ <= 1 and $seen{$_} //= 0 } @_;
	sort {$a <=> $b } keys %seen;
}

# bounding box
sub bbox {
	my $s = shift;
	return $s->{bbox} if $s->{bbox};
	given(order($s)) {
		when(1) {	
			my @p = @{$s->{pt}};
			$s->{bbox} =
				range($p[0][0], $p[0][0], $p[0][1], $p[0][1]) |
				range($p[1][0], $p[1][0], $p[1][1], $p[1][1]);
		}
		default {
			my $pt = $s->{pt}[0];
			my $poly = $s->poly_form(-$pt);
			my $dev = $poly->derivative;

			my @roots = valid_param(1,
					roots([map {$_->[0]} @$dev]),
					roots([map {$_->[1]} @$dev])
				);

			$s->{bbox} = range($pt->[0], $pt->[0], $pt->[1], $pt->[1]);
			for (@roots) {
				$s->{bbox}->include(
					$poly->poly_eval($_) + $pt
				);
			}
		}
	}
	$s->{bbox}
}

sub extrema {
	my $seg = shift;
	my $dv1 = $seg->poly_form(-($seg->{pt}[0]))->derivative;
	valid_param(0, 1,
		poly(map {$_->[0] } @$dv1)->roots,
		poly(map {$_->[1] } @$dv1)->roots
	)
}

# bezier->distance(point)
# calculate distance from point to nearest point on curve by finding
# roots for dot_product(B - x, B')
sub distance {
	my ($seg, $p, $near) = @_;

	my $poly = $seg->poly_form(-$p);
	my $deriv = $poly->derivative;

	my $xs = poly(map { $_->[0] } @$poly) * poly(map {$_->[0] } @$deriv);
	my $ys = poly(map { $_->[1] } @$poly) * poly(map {$_->[1] } @$deriv);
	my @r = grep { !ref $_} (($xs + $ys)->roots);
	return unless @r;
	if (defined $near) {
		return ((sort { abs($a-$near) <=> abs($b-$near) } @r)[0])
	} else {
		return ((sort { $a->[0] <=> $b->[0] }
			map { $seg->param_to_point($_)->dist2($p) }
			@r)[0]->[1])
	}
}

# takes n parameters, and return n-1 subdivided curves
# caller is responsible of ordering and range of the parameters
sub subdivide {
	my $s = shift;
	my $p = poly_form($s);
	my (@pts, @ret) = map { $_->[1]//=$s->param_to_point($_->[0]) } @_;
	my @parent = @{$s->{parent}};
	for my $i (0 .. ($#_ - 1)) {
		my $t0 = $_[$i][0];
		my $dt = $_[$i + 1][0] - $t0;
		my $seg;
		given($#$p) {
		when(1) { $seg = path_segment($pts[$i], $pts[$i+1]) }
		when(2) {
			$seg = from_poly([
				$pts[$i],
				($p->[1] + 2 * $t0 * $p->[2]) * $dt,
				$p->[2] * ($dt * $dt)])
		}
		when(3) {
			my @poly = (
				$pts[$i],
				($p->[1] + 2 * $t0 * $p->[2] + 3 * $t0 * $t0 * $p->[3]) * $dt,
				($p->[2] + 3 * $t0 * $p->[3]) * ($dt * $dt),
				$p->[3] * ($dt * $dt * $dt)
			);

			while ($poly[-1]->len2 < 1e-10) {
				pop @poly;
			}

			next if @poly < 2;
			$seg = from_poly(\@poly);
		}}

		$seg->{pt}[0] = $pts[$i];
		$seg->{pt}[-1] = $pts[$i+1];
		$seg->{parent} = [ $parent[0],
				$parent[1] + ($parent[2] - $parent[1]) * $t0,
				$parent[1] + ($parent[2] - $parent[1]) * ($t0 + $dt)];
		push @ret, $seg;
	}
	#print STDERR "div $s->{uid}: @{[ map {$_->{uid}} @ret ]}";
	wantarray ? @ret : $ret[-1];
}

sub split {
	my $seg = shift;
	return $seg unless $seg->{sect};
	my @split = ([0, $seg->{pt}[0]]);
	for (sort { $a->[0] <=> $b->[0] } @{$seg->{sect}}) {
		if ((my $d = $_->[0] - $split[-1][0]) < 1e-5) {
			printf STDERR "WARN: Very small segment split value %g for $seg\n", $d if $^W && $d;
			next;
		}
		next if $_->[0] >= 1;
		push @split, $_;
	}

	if ($split[-1][0] > 1-1e-8) {
		$split[-1] = [1, $seg->{pt}[-1]]
	} else {
		push @split, [1, $seg->{pt}[-1]]
	}
	if (@split == 2) {
		$seg->{sect} = [];
		return ($seg);
	}
	$seg->subdivide(@split)
}

sub section {
	my ($seg, $range) = @_;
	if (ref $range) {
		return subdivide($seg, [$range->[0]], [$range->[1]])
	} else {
		return $seg->param_to_point($range)
	}
}

# if cubic curve crosses self, return two parameters n1 and n2
# that are at the same point; n1 < n2
sub self_intersect {
	my $seg = shift;
	# if not cubic, return null
	return unless @{$seg->{pt}} >= 4;

	my $p = poly_form($seg);
	my ($a12, $a13, $a23) = ($p->[1]^$p->[2], $p->[1]^$p->[3], $p->[2]^$p->[3]);
	return unless $a23;

	my $plus = - $a13 / $a23;
	# if intersections are out of [0, 1] range
	return unless $plus > 0;

	my $quad = $a12 / $a23;
	my $minus2 = 4 * $quad - 3 * $plus * $plus;
	return unless $minus2 > 0;

	$minus2 = sqrt($minus2);
	my @r = (($plus - $minus2)/2, ($plus + $minus2) / 2);
	return if ($r[0] < 0 || $r[1] > 1);
	my $pt = $seg->param_to_point($r[0]);
	for (@r) {
		next if $_ > 1 || $_ < 0;
		push @{$seg->{sect}}, [$_, $pt]
	}
}

sub sect_11 {
	my ($this, $that, $nolimit) = @_;
	my ($a, $b) = ($this->poly_form, $that->poly_form);
	my ($d, $x) = ($a->[0] - $b->[0], $b->[1] ^ $a->[1]);
	return unless abs($x) > 1e-10;

	my @ret = ($d^$b->[1] / $x, $d^$a->[1] / $x);
	my $tol = 1e-8;

	for (@ret) {
		return unless $_ >= -$tol && $_ <= 1+$tol
	}
	\@ret
}

sub sect_13 {
	my ($this, $that) = @_;

	my $off = -($that->{pt}[0]);
	my $p1 = $this->poly_form($off);
	my ($a0, $b0, $a1, $b1) = map { @$_ } @$p1;
	my ($c0, $d0, $c1, $d1, $c2, $d2, $c3, $d3) = 
		map { @$_ } @{$that->poly_form($off)};

	$c3 //= 0;
	$d3 //= 0;
	my @r;
	for my $s (valid_param roots([
		$a1 * $d0 - $b1 * $c0 - $a1 * $b0 + $a0 * $b1,
		$a1 * $d1 - $b1 * $c1,
		$a1 * $d2 - $b1 * $c2,
		$a1 * $d3 - $b1 * $c3,]
	))
	{
		my $x = $that->param_to_point($s);
		my $t = (($x - $p1->[0] + $off) & $p1->[1]) / ($p1->[1] & $p1->[1]);
		if ($t >= 0 && $t <= 1) {
			push @r, [$t, $s];
		}
	}
	return @r;
}

sub sect_22 {
	my ($this, $that) = @_;
	my $off = -($this->{pt}[0]);

	my @pt = @{$that->poly_form($off)};
	my ($c0, $c1, $c2) = map {$_->[0]} @pt;
	my ($d0, $d1, $d2) = map {$_->[1]} @pt;

	# this new @pt will be used later
	@pt = @{$this->poly_form($off)};
	my ($a0, $a1, $a2) = map {-$_->[0]} @pt;
	my ($b0, $b1, $b2) = map {-$_->[1]} @pt;

	$c0 += $a0;
	$d0 += $b0;

	my $ab12   = $a1*$b2 - $a2*$b1;
	my $adbc20 = $a2*$d0 - $b2*$c0;
	my $bcad21 = $b2*$c1 - $a2*$d1;
	my $bcad12 = $b1*$c2 - $a1*$d2;
	my $bcad22 = $b2*$c2 - $a2*$d2;

	my @coef = (
		($b1*$c0-$a1*$d0)*$ab12 - $adbc20 * $adbc20,
		($b1*$c1-$a1*$d1)*$ab12 + 2*$bcad21*$adbc20,
		-$bcad21*$bcad21 + 2*$adbc20*$bcad22 + $bcad12*$ab12,
		-2*$bcad21*$bcad22,
		-$bcad22 * $bcad22,
	);

	my @ret;
	for my $s (valid_param roots(\@coef)) {
		my $p = $that->param_to_point($s) + $off;
		my $c = $pt[1]^$pt[2] or next;
		my $t = -(($pt[0] - $p)^$pt[2] / $c);
		next if $t > 1 || $t < 0;
		push @ret, [$t, $s];
	}
	@ret;
}

sub sect_33 {
	my ($this, $that) = @_;
	my $off = -($this->{pt}[0]);

	my @pt = @{$that->poly_form($off)};
	my ($c0, $c1, $c2, $c3) = map {$_->[0]} @pt;
	my ($d0, $d1, $d2, $d3) = map {$_->[1]} @pt;

	# this new @pt will be used later
	@pt = @{$this->poly_form($off)};
	my ($a0, $a1, $a2, $a3) = map {$_->[0]} @pt;
	my ($b0, $b1, $b2, $b3) = map {$_->[1]} @pt;

	$a3 //= 0;
	$b3 //= 0;
	$c3 //= 0;
	$d3 //= 0;

	my $ac = $a0 - $c0;
	my $bd = $b0 - $d0;
	my $m33 = $b3*$c3-$a3*$d3;
	my $m32 = $b3*$c2-$a3*$d2;
	my $m23 = $b2*$c3-$a2*$d3;
	my @coef = (
		$a3*$a3*$a3*$bd*$bd*$bd-$a3*$a3*($a2*$b1*$bd*$bd+$a1*$bd*(
		-$b1*$b1+2*$b2*$bd)+$ac*($b1*$b1*$b1+3*$b3*$bd*$bd-3*$b1*$b2
		*$bd))-$b3*($a2*$a2*$a2*$bd*$bd+$b3*($a1*$a1*$b1*$ac-$a1*$b2
		*$ac*$ac+$b3*$ac*$ac*$ac-$a1*$a1*$a1*$bd)+$a2*$a2*($a1*$b1*(
		-$bd)+$ac*($b1*$b1-2*$b0*$b2+2*$b2*$d0))+$a2*(($b2*$b2-2*$b1
		*$b3)*$ac*$ac+$a1*$a1*$b2*$bd-$a1*$ac*($b1*$b2-3*$b0*$b3+3
		*$b3*$d0)))+$a3*($a1*$ac*(-$b1*$b2*$b2+2*$b1*$b1*$b3+$b2*$b3
		*$bd)+$ac*$ac*($b2*$b2*$b2-3*$b1*$b2*$b3+3*$b3*$b3*$bd)+$a1
		*$a1*($b2*$b2-2*$b1*$b3)*$bd+$a2*$a2*$b2*$bd*$bd+$a2*($a1*(
		-$b1*$b2+3*$b3*$bd)*$bd+$ac*($b1*$b1*$b2-2*$b2*$b2*$bd-$b1
		*$b3*$bd))),

		-3*$a3*$a3*$a3*$bd*$bd*$d1+$a3*$a3*($b1*$b1*$b1*$c1-$a1*$b1
		*$b1*$d1+$bd*(3*$b0*$b3*$c1-3*$b3*$c1*$d0+4*$a1*$b2*$d1+6
		*$a0*$b3*$d1-6*$b3*$c0*$d1)+$b1*(-3*$b0*$b2*$c1+3*$b2*$c1
		*$d0+2*$a2*$b0*$d1-3*$a0*$b2*$d1+3*$b2*$c0*$d1-2*$a2*$d0
		*$d1))+$a3*(-($a1*$a1*($b2*$b2-2*$b1*$b3)*$d1)-2*$a2*$a2*$b2
		*$bd*$d1-$ac*(2*$b2*$b2*$b2*$c1-6*$b1*$b2*$b3*$c1+3*$b3*$b3
		*(2*$b0*$c1-2*$c1*$d0+$a0*$d1-$c0*$d1))+$a1*($b1*$b2*$b2*$c1
		-2*$b1*$b1*$b3*$c1+$b2*$b3*(-$b0*$c1+$c1*$d0-$a0*$d1+$c0
		*$d1))+$a2*(-$b1*$b1*$b2*$c1+2*(3*$a1*$b3*$d0*$d1+$b0*($b2
		*$b2*$c1-3*$a1*$b3*$d1)-$b2*$b2*($c1*$d0-$a0*$d1+$c0*$d1))
		+$b1*($b0*$b3*$c1+$a1*$b2*$d1-$b3*($c1*$d0-$a0*$d1+$c0
		*$d1))))+$b3*(2*$a2*$a2*$a2*$bd*$d1+$b3*($a1*$a1*$b1*$c1+3
		*$b3*$ac*$ac*$c1-2*$a1*$b2*$ac*$c1-$a1*$a1*$a1*$d1)+$a2*$a2
		*($b1*$b1*$c1-$a1*$b1*$d1+2*$b2*(-$b0*$c1+$c1*$d0-$a0*$d1
		+$c0*$d1))+$a2*(2*($b2*$b2-2*$b1*$b3)*$ac*$c1+$a1*$a1*$b2
		*$d1-$a1*($b1*$b2*$c1+3*$b3*(-$b0*$c1+$c1*$d0-$a0*$d1+$c0
		*$d1)))),

		-3*$a3*$a3*$a3*$bd*(-$d1*$d1+$b0*$d2-$d0*$d2)+$a3*$a3*($b1
		*$b1*$b1*$c2+3*$b0*$b0*$b3*$c2+3*$b3*$c2*$d0*$d0+6*$b3*$c1
		*$d0*$d1-2*$a1*$b2*$d1*$d1-3*$a0*$b3*$d1*$d1+3*$b3*$c0*$d1
		*$d1-$a1*$b1*$b1*$d2-4*$a1*$b2*$d0*$d2-6*$a0*$b3*$d0*$d2+6
		*$b3*$c0*$d0*$d2+$b0*(4*$a1*$b2*$d2-6*$b3*($c2*$d0+$c1*$d1
		-$a0*$d2+$c0*$d2))+$b1*($b0*(-3*$b2*$c2+2*$a2*$d2)+3*$b2
		*($c2*$d0+$c1*$d1-$a0*$d2+$c0*$d2)-$a2*($d1*$d1+2*$d0*$d2)))
		+$a3*($b2*$b2*$b2*($c1*$c1-2*$ac*$c2)+$b2*$b2*($a1*($b1*$c2
		-$a1*$d2)+2*$a2*($b0*$c2-$c2*$d0-$c1*$d1+$a0*$d2-$c0*$d2))
		+$b2*(-$a2*$b1*$b1*$c2+$b1*(-3*$b3*($c1*$c1-2*$a0*$c2+2*$c0
		*$c2)+$a1*$a2*$d2)+$a1*$b3*(-$b0*$c2+$c2*$d0+$c1*$d1-$a0*$d2
		+$c0*$d2)+$a2*$a2*($d1*$d1-2*$b0*$d2+2*$d0*$d2))+$b3*(-2*$a1
		*$b1*$b1*$c2-3*$b3*$c1*$c1*$d0-$a2*$b1*$c2*$d0+6*$a0*$b3*$c2
		*$d0-6*$b3*$c0*$c2*$d0-$a2*$b1*$c1*$d1+6*$a0*$b3*$c1*$d1-6
		*$b3*$c0*$c1*$d1+3*$a1*$a2*$d1*$d1+2*$a1*$a1*$b1*$d2+$a0*$a2
		*$b1*$d2-3*$a0*$a0*$b3*$d2-$a2*$b1*$c0*$d2+6*$a0*$b3*$c0*$d2
		-3*$b3*$c0*$c0*$d2+6*$a1*$a2*$d0*$d2+$b0*(3*$b3*($c1*$c1-2
		*$a0*$c2+2*$c0*$c2)+$a2*($b1*$c2-6*$a1*$d2))))+$b3*($b3*($a1
		*$a1*$b1*$c2+3*$b3*$ac*(-$c1*$c1+$a0*$c2-$c0*$c2)+$a1*$b2
		*($c1*$c1-2*$a0*$c2+2*$c0*$c2)-$a1*$a1*$a1*$d2)-$a2*$a2*$a2
		*($d1*$d1-2*$bd*$d2)+$a2*$a2*($b1*$b1*$c2-$a1*$b1*$d2+2*$b2
		*(-$b0*$c2+$c2*$d0+$c1*$d1-$a0*$d2+$c0*$d2))-$a2*($b2*$b2
		*($c1*$c1-2*$ac*$c2)+$a1*$b2*($b1*$c2-$a1*$d2)+$b3*(-2*$b1
		*($c1*$c1-2*$a0*$c2+2*$c0*$c2)+3*$a1*(-$b0*$c2+$c2*$d0+$c1
		*$d1-$a0*$d2+$c0*$d2)))),

		$b3*$b3*$b3*($c1*$c1*$c1-6*$ac*$c1*$c2+3*$ac*$ac*$c3)+$b3
		*$b3*($a1*(2*$b2*($c1*$c2-$ac*$c3)+$a1*($b1*$c3-$a1*$d3))-3
		*$a3*(-2*$b0*($c1*$c2-$ac*$c3)+$c1*$c1*$d1+2*$c1*($c2*$d0
		-$a0*$d2+$c0*$d2)+$ac*(-2*$c3*$d0-2*$c2*$d1+$a0*$d3-$c0
		*$d3))+$a2*(4*$b1*($c1*$c2-$ac*$c3)+3*$a1*($b0*$c3-$c3*$d0
		-$c2*$d1-$c1*$d2+$a0*$d3-$c0*$d3)))+$b3*(3*$a3*$a3*($b0*$b0
		*$c3+$c3*$d0*$d0+2*$c2*$d0*$d1+$c1*$d1*$d1+2*$c1*$d0*$d2-2
		*$a0*$d1*$d2+2*$c0*$d1*$d2-2*$a0*$d0*$d3+2*$c0*$d0*$d3-2*$b0
		*($c3*$d0+$c2*$d1+$c1*$d2-$a0*$d3+$c0*$d3))+$a3*(-2*$a1*$b1
		*$b1*$c3-$b1*(6*$b2*($c1*$c2-$ac*$c3)-2*$a1*$a1*$d3+$a2*(
		-$b0*$c3+$c3*$d0+$c2*$d1+$c1*$d2-$a0*$d3+$c0*$d3))+$a1*(
		-($b0*($b2*$c3+6*$a2*$d3))+$b2*($c3*$d0+$c2*$d1+$c1*$d2-$a0
		*$d3+$c0*$d3)+6*$a2*($d1*$d2+$d0*$d3)))+$a2*(-2*$b2*$b2*($c1
		*$c2-$ac*$c3)+$b2*($a1*(-$b1*$c3+$a1*$d3)+2*$a2*(-$b0*$c3
		+$c3*$d0+$c2*$d1+$c1*$d2-$a0*$d3+$c0*$d3))+$a2*($b1*$b1*$c3
		-$a1*$b1*$d3-2*$a2*($d1*$d2-$b0*$d3+$d0*$d3))))-$a3*(-2*$b2
		*$b2*$b2*($c1*$c2-$ac*$c3)+$b2*$b2*($a1*(-$b1*$c3+$a1*$d3)+2
		*$a2*(-$b0*$c3+$c3*$d0+$c2*$d1+$c1*$d2-$a0*$d3+$c0*$d3))+$a3
		*(-$b1*$b1*$b1*$c3+$a1*$b1*$b1*$d3+$a3*($d1*$d1*$d1-6*$bd
		*$d1*$d2+3*$bd*$bd*$d3)+2*$a2*$b1*($d1*$d2-$bd*$d3))+$b2
		*($a3*($b0*(3*$b1*$c3-4*$a1*$d3)-3*$b1*($c3*$d0+$c2*$d1+$c1
		*$d2-$a0*$d3+$c0*$d3)+4*$a1*($d1*$d2+$d0*$d3))+$a2*($b1*$b1
		*$c3-$a1*$b1*$d3-2*$a2*($d1*$d2-$bd*$d3)))),

		3*$b3*$b3*$b3*($c1*$c1*$c2-$ac*$c2*$c2-2*$ac*$c1*$c3)+$b3
		*$b3*($a1*$b2*($c2*$c2+2*$c1*$c3)+3*$a3*($b0*($c2*$c2+2*$c1
		*$c3)-$c2*$c2*$d0-2*$c1*$c3*$d0-2*$c1*$c2*$d1+2*$a0*$c3*$d1
		-2*$c0*$c3*$d1-$c1*$c1*$d2+2*$a0*$c2*$d2-2*$c0*$c2*$d2+2*$a0
		*$c1*$d3-2*$c0*$c1*$d3)+$a2*(2*$b1*($c2*$c2+2*$c1*$c3)-3*$a1
		*($c3*$d1+$c2*$d2+$c1*$d3)))+$b3*(3*$a3*$a3*(2*$c3*$d0*$d1
		+$c2*$d1*$d1+2*$c2*$d0*$d2+2*$c1*$d1*$d2-$a0*$d2*$d2+$c0*$d2
		*$d2+2*$c1*$d0*$d3-2*$a0*$d1*$d3+2*$c0*$d1*$d3-2*$b0*($c3
		*$d1+$c2*$d2+$c1*$d3))-$a2*($b2*$b2*($c2*$c2+2*$c1*$c3)-2
		*$a2*$b2*($c3*$d1+$c2*$d2+$c1*$d3)+$a2*$a2*($d2*$d2+2*$d1
		*$d3))+$a3*(-($b1*(3*$b2*$c2*$c2+6*$b2*$c1*$c3+$a2*$c3*$d1
		+$a2*$c2*$d2+$a2*$c1*$d3))+$a1*($b2*$c3*$d1+$b2*$c2*$d2+3
		*$a2*$d2*$d2+$b2*$c1*$d3+6*$a2*$d1*$d3)))+$a3*($b2*$b2*$b2
		*($c2*$c2+2*$c1*$c3)-2*$a2*$b2*$b2*($c3*$d1+$c2*$d2+$c1*$d3)
		-$a3*($a2*$b1*($d2*$d2+2*$d1*$d3)+3*$a3*($d1*$d1*$d2-$bd*$d2
		*$d2-2*$bd*$d1*$d3))+$b2*($a2*$a2*($d2*$d2+2*$d1*$d3)+$a3*(3
		*$b1*($c3*$d1+$c2*$d2+$c1*$d3)-2*$a1*($d2*$d2+2*$d1*$d3)))),

		3*$b3*$b3*$b3*($c1*$c2*$c2+$c1*$c1*$c3-2*$ac*$c2*$c3)+$b3*(
		-2*$a2*(-$b2*$c2+$a2*$d2)*(-$b2*$c3+$a2*$d3)+3*$a3*$a3*($c1
		*$d2*$d2+$c3*($d1*$d1-2*$bd*$d2)+2*$c1*$d1*$d3-2*$a0*$d2*$d3
		+2*$c0*$d2*$d3+2*$c2*($d1*$d2-$b0*$d3+$d0*$d3))+$a3*(-($b1
		*(6*$b2*$c2*$c3+$a2*$c3*$d2+$a2*$c2*$d3))+$a1*($b2*$c3*$d2
		+$b2*$c2*$d3+6*$a2*$d2*$d3)))-$a3*(-2*$b2*$b2*$b2*$c2*$c3+2
		*$a2*$b2*$b2*($c3*$d2+$c2*$d3)-$b2*(2*$a2*$a2*$d2*$d3+$a3*(3
		*$b1*$c3*$d2+3*$b1*$c2*$d3-4*$a1*$d2*$d3))+$a3*(2*$a2*$b1
		*$d2*$d3+3*$a3*($d1*$d2*$d2+$d1*$d1*$d3-2*$bd*$d2*$d3)))+$b3
		*$b3*(4*$a2*$b1*$c2*$c3+2*$a1*$b2*$c2*$c3-3*$a1*$a2*$c3*$d2
		-3*$a1*$a2*$c2*$d3+$a3*(6*$b0*$c2*$c3-3*($c2*$c2*$d1+2*$c1
		*$c3*$d1-2*$ac*$c3*$d2+$c1*$c1*$d3+2*$c2*($c3*$d0+$c1*$d2
		-$a0*$d3+$c0*$d3)))),

		$b3*$b3*$b3*($c2*$c2*$c2+6*$c1*$c2*$c3-3*$ac*$c3*$c3)+$b3
		*$b3*($c3*(2*$a2*$b1*$c3+$a1*$b2*$c3-3*$a1*$a2*$d3)+3*$a3
		*($b0*$c3*$c3-$c3*$c3*$d0-2*$c3*($c2*$d1+$c1*$d2-$a0*$d3+$c0
		*$d3)-$c2*($c2*$d2+2*$c1*$d3)))+$b3*(-$a2*$m23*$m23+$a3*(
		-($b1*$c3*(3*$b2*$c3+$a2*$d3))+$a1*$d3*($b2*$c3+3*$a2*$d3))
		+3*$a3*$a3*($d3*(2*$c1*$d2-$a0*$d3+$c0*$d3)+2*$c3*($d1*$d2
		-$bd*$d3)+$c2*($d2*$d2+2*$d1*$d3)))-$a3*(-$b2*$b2*$b2*$c3
		*$c3+2*$a2*$b2*$b2*$c3*$d3+$b2*$d3*(-3*$a3*$b1*$c3-$a2*$a2
		*$d3+2*$a1*$a3*$d3)+$a3*($a2*$b1*$d3*$d3+$a3*($d2*$d2*$d2+6
		*$d1*$d2*$d3-3*$bd*$d3*$d3))),

		3*$m33*($b3*$b3*($c2*$c2+$c1*$c3)-$a3*$b3*($c3*$d1+2*$c2*$d2
		+$c1*$d3)+$a3*$a3*($d2*$d2+$d1*$d3)),

		3*$m32*$m33*$m33,

		$m33*$m33*$m33,

	);

	return unless grep { abs($_) > 1e-9 } @coef;

	my @crosses = ();
	for my $s (valid_param roots(\@coef)) {
		my $x = $that->param_to_point($s) + $off;
		my ($a31, $a32, $a1x, $a2x, $a3x) = (
			$pt[3]^$pt[1],
			$pt[3]^$pt[2],
			$pt[1]^($pt[0] - $x),
			$pt[2]^($pt[0] - $x),
			$pt[3]^($pt[0] - $x),
		);

		my ($num, $denom) = (-($a1x * $a31 - $a2x * $a3x), $a32 * $a1x - $a3x * $a3x);
		next if (abs($denom) < 1e-10);
		my $t = $num / $denom;
		if (0<= $t && $t <= 1) {
			push @crosses, [$t, $s];
		}
	}
	return @crosses;
}

my %sects = (
	'11'	=> [\&sect_11, 0],
	'12'	=> [\&sect_13, 0],
	'21'	=> [\&sect_13, 1],
	'22'	=> [\&sect_22, 0],
	'13'	=> [\&sect_13, 0],
	'31'	=> [\&sect_13, 1],
	'32'	=> [\&sect_33, 0],
	'23'	=> [\&sect_33, 1],
	'33'	=> [\&sect_33, 0],
);

sub intersect {
	my ($this, $that, $must_intersect) = @_;
	my $c = $this->bbox & $that->bbox;
	return unless defined $c;

	if (exists $that->{kids}) {
		for (@{$that->{kids}}) {
			$this->intersect($_)
		}
		return;
	}
	my $key = order($this).order($that);
	if ($key eq '33') {
		return if $c->[0]==$c->[1]||$c->[2]==$c->[3];
	}

	my @ret;
	unless($must_intersect) {
		@ret = overlap($this, $that)
	}
	if (!@ret) {
		my ($func, $swap) = @{$sects{$key}};
		@ret = ($swap	? $func->($that, $this)
				: $func->($this, $that)
			);

		if (!@ret) {
			return unless $must_intersect;
			die "Failed to intersect $this, $that";
		};

		@ret = map { fix_sect($this, $that, @$_) }
			($swap	? map { [$_->[1], $_->[0]] } @ret
				: @ret);
	}
	if (!$must_intersect) {
		for (@ret) {
			push @{$this->{sect}}, $_->[0]
				if ($_->[0][0] > 0 && $_->[0][0] < 1);
			push @{$that->{sect}}, $_->[1]
				if ($_->[1][0] > 0 && $_->[1][0] < 1);
		}
		return
	}
	for (@ret) {
		next unless 	$_->[0][0] > 0 &&
				$_->[0][0] < 1 &&
				$_->[1][0] > 0 &&
				$_->[1][0] < 1;

		push @{$this->{sect}}, [$_->[0][0], $_->[1][0], $that];
	}
}

sub fix_sect {
	my ($this, $that, $t1, $t2) = @_;
	my ($ot1, $ot2, $p1, $p2) = ($t1, $t2,
		$this->param_to_point($t1), $that->param_to_point($t2));

	my $dist = $p1->dist2($p2);
	my $odist = $dist;

	my $tries = 0;
	while ($dist > 1e-10) {
		$p1 = $this->param_to_point($t1);
		$t2 = ($t2 + $that->distance($p1, $t2)) / 2;

		$p2 = $that->param_to_point($t2);
		$t1 = ($t1 + $this->distance($p2, $t1)) / 2;

		$dist = $p1->dist2($p2);
		if (++$tries > 20) {
			carp join("\n", @_), "\nroot still not accurate after too many tries; dist: $dist, ($p1), ($p2)" if $dist > 1e-6;
			last;
		}
	}
	if ($dist > $odist) {
		$p1 = $this->param_to_point($ot1);
		print STDERR "using old value";
		return [[$ot1, $p1], [$ot2, $p2]]
	}
	[ [$t1, $p1], [$t2, $p1] ];
}

# finds out if curve a(t) can be express via poly_form of curve b(s) with param
# transformation t = s0 + alpha * s
sub overlap3 {
	my ($a, $b) = map { poly_form $_ } @_;

	my $alpha = vdiv($b->[3], $a->[3]);
	return unless $alpha && (abs($alpha) > 1e-12);
	$alpha = ($alpha > 0) ? $alpha**(1/3) : -(-$alpha)**(1/3);
	my $s0 = vdiv($b->[2] / ($alpha*$alpha) - $a->[2], $a->[3]);
	return unless	defined $s0;
	$s0 /= 3;
	return unless
		veq($b->[1], $alpha * ($a->[1] + 2*$s0*$a->[2] + 3*$s0*$s0*$a->[3])) and
		veq($b->[0], $a->poly_eval($s0));

	return ($alpha, $s0);
}

sub overlap2 {
	my ($a, $b) = map { poly_form $_} @_;

	my $alpha = vdiv($b->[2], $a->[2]);
	return unless $alpha && $alpha > 1e-12;

	$alpha = sqrt $alpha;
	my $s0 = vdiv($b->[1] / $alpha - $a->[1], $a->[2]);
	if (! defined $s0) {
		$alpha = - $alpha;
		$s0 = vdiv($b->[1] / $alpha - $a->[1], $a->[2]);
	}
	return unless	defined $s0;

	$s0 /= 2;
	return unless	veq($b->[0], $a->poly_eval($s0));

	return ($alpha, $s0);
}

sub overlap1 {
	my ($a, $b) = map { poly_form $_ } @_;
	my $alpha = vdiv($b->[1], $a->[1]);

	return unless $alpha;
	my $s0 = vdiv($b->[0] - $a->[0], $a->[1]);
	return defined($s0) ? ($alpha, $s0) : ();
}

sub overlap {
	my ($this, $that) = @_;
	my ($ra, $rb, $alpha, $s0);

	if ($this->{parent}[0] == $that->{parent}[0]) {
		my @r1 = @{$this->{parent}}[1,2];
		my @r2 = @{$that->{parent}}[1,2];
		my @r = (max(min(@r1), min(@r2)), min(max(@r1), max(@r2)));

		if ($r[0] >= $r[1]) { return }

		if (($r2[1]-$r2[0]) * ($r1[1]-$r1[0]) < 0) {
			@r = ($r[1], $r[0]);
		}

		$s0 = ($r[0] - $r1[0]) / ($r1[1]-$r1[0]);
		$alpha = ($r[1] - $r[0]) / ($r1[1] - $r1[0]);
		print STDERR"same @r1 @r2 $s0 $alpha";
	} else {
		return unless order($this) == order($that) &&
			$this->bbox & $that->bbox;

		my $func;
		given (order($this)) {
			when(1) { $func = \&overlap1 }
			when(2) { $func = \&overlap2 }
			when(3) { $func = \&overlap3 }
		}
		($alpha, $s0) = $func->($this, $that) or return;
	}

	$ra = [$s0, $s0 + $alpha, $this];
	$rb = $alpha < 0 
			? [(1-$s0)/$alpha, -$s0/$alpha, $that]
			: [-$s0/$alpha, (1-$s0)/$alpha, $that];
	my ($p1, $p2);
	for ($ra, $rb) {
		if ($_->[0] < 0) {
			$_->[0] = 0;
			$p1 = $_->[2]{pt}[0];
		}
		if ($_->[0] > 1) {
			$_->[0] = 1;
			$p1 = $_->[2]{pt}[-1];
		}
		if ($_->[1] < 0) {
			$_->[1] = 0;
			$p2 = $_->[2]{pt}[0];
		}
		if ($_->[1] > 1) {
			$_->[1] = 1;
			$p2 = $_->[2]{pt}[-1];
		}
	}
	return if $ra->[0] == $ra->[1] or $rb->[0] == $rb->[1];

	$p1 //= $this->param_to_point($ra->[0]);
	$p2 //= $this->param_to_point($ra->[1]);
	(	[[$ra->[0], $p1], [$rb->[0], $p1]],
		[[$ra->[1], $p2], [$rb->[1], $p2]]
	);
}

sub derivative {
	my $seg = shift;
	my $poly = $seg->poly_form(-($seg->{pt}[0]))->derivative;
	map { $poly->poly_eval($_) } @_
}

# segment->(parameter, limit_from_left)
sub tangent {
	my $p = shift;
	my ($s, $t) = ($p->poly_form->derivative, shift);
	my $left = shift;

	#short lines
	my $x = $s->poly_eval($t);
	if ($x->len2 > 1e-6 || $p->order == 1) {
		return $x->normalize;
	}
	$s = $s->derivative;
	$x = $s->poly_eval($t)
		or confess "Bad tangent: $p at $t, from left: $left";

	if ($x->len2 > 1e-6) {
		return ($left ? -$x : $x)->normalize;
	}
	$s = $s->derivative;
	return vector(0,0) if (!@$s);
	return $s->poly_eval($t)->normalize
}

sub flat {
	my $s = shift;
	return $s unless order($s) > 1;
	my @p = @{$s->poly_form};
	for my $i (1 .. ($#p-1)) {
		return $s if abs($p[$i]^$p[$i + 1]) > 1e-6
	}
	my ($p1, $p2) = ($s->{pt}[0], $s->{pt}[-1]);
	if ($p1->dist2($p2) > 1e-7) {
		return path_segment($s->{pt}[0], $s->{pt}[-1]);
	}
}

sub transform {
	my ($ps, $m) = @_;
	path_segment(map { $_->transform($m) } @{$ps->{pt}})
}

sub to_cubic {
	my @pt = @{shift->{pt}};
	@pt == 3 or return;
	path_segment(	$pt[0], ($pt[1]*2 + $pt[0])/3,
			($pt[1]*2 + $pt[2])/3, $pt[2]	);
}

1;

